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Resumen 


De manera tradicional, los sistemas de información geográfica disponen de 
información sobre el territorio en dos formatos: formato ráster y formato vectorial, 
y cada uno de estos formatos es tratado con sus algoritmos independientes. Por 
consiguiente, cualquier tipo de estudio hidráulico realizado sobre el territorio es 
una información susceptible de ser clasificada en estos dos formatos y también 
dispone de sus propios algoritmos. Por otro lado, LiDAR es la tecnología más precisa 
para generar Modelos Digitales de Elevaciones (MDE) para grandes áreas, con una 
resolución espacial de 1 0 2 m, y con una gran precisión en altura. Ello confiere un 
carácter ráster a la información cartográfica, cuya unidad básica de información 
es la celda. Por consiguiente, este tipo de material resulta muy adecuado para la 
generación de mallas de cuadriláteros, estructuradas para los modelos numéricos 
de flujo en lámina libre 2D. El objetivo principal es presentar una adaptación del 
esquema clásico explícito de primer orden de volúmenes finitos a una malla de 
Volúmenes Finitos Cuadrangulares (VFC), que hace equivaler —directamente y sin 
interpolaciones— cada celda del MDE (o ráster), con el correspondiente VFC. Por 
consiguiente, las dimensiones del VFC serán las de la Celda del MDE de base. Ello 
se presenta aquí como el enfoque ráster al problema del problema hidrodinámico del flujo 
en lámina libre en 2D. En el enfoque ráster, el hecho de que a cada instante de tiempo 
de integración no sea necesaria la consulta de la base de datos de la topología de 
la malla hace que el proceso numérico mejore de forma enorme su eficiencia de 
cálculo. Por consiguiente, se pueden abordar problemas con mayor dimensión. 
También es objetivo comparar ambos enfoques a través de la resolución de un 
ejemplo ilustrativo. 


Palabras clave: modelización bidimensional, paralelización, aguas someras, 
volúmenes finitos. 


Introducción y objetivos 


La información del territorio con la que los 
sistemas de información geográfica (SIG) tra- 
bajan viene dada en dos formatos: el formato 
vectorial y el formato ráster. En las plataformas 
de SIG existen algoritmos independientes para 
el tratamiento de la información, dependiendo 
del formato original de los datos. Por lo tanto, 
cualquier información procedente de un modelo 


hidráulico de flujo de agua en lámina libre 
circulando por el territorio puede ser clasificada 
en estos dos formatos. Por consiguiente, los 
algoritmos de simulación también pueden ser 
independientes, con lo que son susceptibles 
de ser clasificados de dos maneras: enfoque 
vectorial y enfoque ráster. 

La mayoría de los modelos numéricos de 
simulación del flujo en lámina libre en 2D se 
basan en unas representaciones geométricas del 
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territorio, denominadas mallas. Dichas mallas 
están compuestas de unidades simples, llama- 
das celdas o volúmenes finitos, que son polí- 
gonos de tres o cuatro lados. La implementación 
de la malla se hace considerando que en el 
interior de sus celdas se mantienen constantes 
las propiedades geométricas y las hidráulicas. 
Por lo tanto, en el contexto de los SIG, las mallas 
son cubiertas vectoriales de polígonos con sus 
correspondientes propiedades. En este trabajo 
se presenta una herramienta de simulación del 
flujo basada en un tipo de celdas denominadas 
volúmenes finitos cuadrangulares (VFC). 

En el contexto de los modelos numéricos, 
cuando se requiere más precisión en los 
resultados o mayor cantidad de información 
en determinadas zonas del ámbito de estudio, 
es necesario crear mallas de volúmenes finitos 
más densas en dichas zonas y menos densas en 
el resto. Ello da lugar a mallas de volúmenes 
finitos no homogéneos, con lo que el enfoque 
vectorial resulta ser el más adecuado para este 
caso. Pero en determinadas situaciones no es 
necesaria una distribución de la información 
tan irregular, por lo que la geometría se puede 
representar mediante mallas uniformes de 
volúmenes finitos homogéneos. Y si además 
la celda base es cuadrangular, entonces el 
problema pasa a tener un enfoque ráster. 

En este trabajo se define un modelo digital 
de elevaciones (MDE) como una representación 
en formato ráster de la geometría de una 
zona del territorio determinada. Un ráster 
está constituido por un conjunto de celdas 
cuadrangulares. Cada celda representa un área 
cuadrada del plano x/y de una determinada 
longitud de lado —longitud que se denomina 
lado de la celda— y un valor asignado a esta 
área. Cuando dicho valor es la cota media 
de toda la superficie de la celda, el ráster 
se denomina MDE. El enfoque ráster en el 
problema hidrodinámico del flujo en lámina 
libre resulta el adecuado cuando el problema 
tiene grandes dimensiones y la información 
original para la construcción de la malla 
proviene directamente de un MDE. Tal es 
el caso de los estudios de inundabilidad o 


modelos para la clasificación de la rotura de un 
embalse en función del riesgo potencial. 

Hoy día existen técnicas numéricas de 
resolución de las ecuaciones de flujo en 
lámina libre en dos dimensiones de Saint 
“Shallow 
Waters Equations”—, mediante esquemas de 


Venant —también denominadas 
volúmenes finitos para todo tipo de mallas, 
bien sean estructuradas y uniformes, bien 
sean no estructuradas y no homogéneas. 
En este trabajo se presenta una adaptación 
de una de estas técnicas numéricas a mallas 
con celdas 
estructuradas de tipo VFC. Dicho en otras 
palabras, se trata de una adaptación numérica 


cuadrangulares uniformes y 


al enfoque ráster porque se establece una malla 
homogénea y estructurada de cuadriláteros 
que hace equivaler —directamente y sin 
interpolaciones— cada celda del MDE con el 
correspondiente VFC. 

Como siempre pasa en hidroinformática, 
cuanto más detallado es el modelo utilizado, 
mejor es el conocimiento que se obtiene del 
fenómeno, pero más difícil es su construcción 
y mayor su complejidad. Antaño, el uso de 
modelos excesivamente simplificados era 
obligada, porque no existían ordenadores con 
la potencia suficiente para el uso de modelos 
sofisticados. Pero hoy en día, la tecnología 
informática disponible permite el uso de 
modelos cada vez más complejos y acordes con 
la realidad, de manera que el establecimiento 
de tales hipótesis simplificadoras resulta 
injustificable. A lo largo de la historia han ido 
apareciendo hipótesis simplificadoras, que en 


la actualidad no tienen por qué ser establecidas: 


e Régimen uniforme. 


+ Régimen estacionario y suavemente 
variable. 

+ Régimen no estacionario y suavemente 
variable. 

+ Unidimensionalidad. 

e  Bidimensionalidad. 

+ Distribución hidrostática de presiones. 


+ Sin pérdidas por infiltración. 
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e Continuidad espacial del coeficiente de 
rugosidad. 

e Existencia de un flujo mínimo como 
condición inicial. 


Por otro lado, la tecnología LIDAR —del 
inglés “Light Detection and Ranging”— es 
una técnica de teledetección que mide la 
altura del terreno utilizando un escáner 
láser. Hoy en día, LiDAR es la tecnología 
más precisa para generar MDE para grandes 
áreas, con una resolución espacial de 0.5 o 1 
m, y con una precisión mínima de 15 cm en 
altura. Esta tecnología ofrece grandes ventajas 
en la obtención de los MDE, en contraposición 
de los sistemas fotogramétricos: permite 
medir la altura real del terreno a través de la 
vegetación, porque el rayo láser la atraviesa, 
tiene una precisión homogénea para toda la 
información del área de estudio, y reduce los 
costes de producción y los plazos de entrega. 
Por consiguiente, este tipo de material resulta 
muy adecuado para la generación de mallas 
estructuradas de VFC uniformes para los 
modelos numéricos de flujo en lámina libre 
2D. 

Por lo general, los modelos numéricos 
comerciales de flujo 2D en lámina libre capa- 
ces de reproducir el comportamiento del 
flujo circulando sobre el territorio requieren 
la generación de complejas mallas para la 
definición geométrica del problema. Estas 
mallas resultan muy tediosas de ser imple- 
mentadas y comprometen gran parte del 
tiempo del personal altamente especializado, 
porque se trata de herramientas adaptadas al 
enfoque vectorial. Además, al crear la malla de 
volúmenes finitos, se pierde precisión por la 
necesidad de interpolación a partir del MDE. 
En resumen, en este trabajo se presenta un 
esquema numérico de volúmenes finitos con 
las siguientes características: 


e Establece una malla estructurada y uni- 
forme de cuadriláteros, que hace equivaler 


—directamente y sin interpolaciones— 
cada celda del MDE con el correspondiente 
VFC. Las dimensiones del VEC serán las 
de la celda de base del MDE, con lo que se 
conserva la precisión original. 

Utiliza el esquema explícito de primer or- 
den en volúmenes finitos de Godunov con 
el Riemann solver de Roe y con el tratamiento 
del término independiente propuesto por 
Vázquez-Cendón (Bladé y Gómez, 2006). 
El uso de un esquema de mayor precisión 
lleva consigo un incremento del tiempo 
de cálculo y no mejora la calidad de los 
resultados, dada la gran precisión de la 
discretización espacial utilizada. 

Esta discretización con malla estructurada 
y uniforme confiere al estudio hidrodiná- 
mico del flujo en dos dimensiones un 
enfoque ráster, en contraposición a la 
tradicional de malla de volúmenes finitos 
no estructurada y no uniforme, que se 
podría decir que tiene un enfoque vectorial. 
El hecho de que a cada instante de tiempo 
de integración no sea necesaria la consulta 
de la base de datos de la topología de 
la malla (nodo, arco y polígono) por ser 
conocida de antemano —porque se trata 
precisamente de una malla estructurada y 
homogénea— hace que el proceso numérico 
mejore de manera enorme su eficiencia de 
cálculo tanto por requerir menor tiempo 
de cálculo como menor memoria para 
el almacenamiento de la geometría. Por 
consiguiente, se pueden abordar problemas 
de mayor dimensión, dado que se consigue 
una mayor optimización del código en aras 
de la obtención de la máxima velocidad de 
cálculo y el uso de información masiva. 

El enfoque ráster es altamente paralelizable 
si se desea. 

Con el enfoque ráster se reduce la comple- 
jidad del preproceso y del postproceso. 
Por contra, con geometrías complicadas, 
donde se requiere la precisión del enfoque 
vectorial, el enfoque ráster resulta ina- 
decuado. 


Tecnología y 
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e La sencillez de planteamiento del enfo- 
que ráster se mantiene también en el 
establecimiento de las condiciones de 
contorno. Se proponen dos tipos de 

condiciones de contorno por defecto: pared 
material sin flujo para los lados contiguos de 
las celdas tipo NODATA —celdas del MDE 
sin información de cota— y salida libre de 
agua —para el contorno del MDE—. 


Después de esta introducción, el presente 
documento se estructura de la siguiente 
manera: la segunda parte de planteamiento 
matemático, donde se describe la discretización 
numérica en volúmenes finitos aplicados 
a las ecuaciones de Saint-Venant y cómo 
las ecuaciones algebraicas resultantes se 
simplifican cuando se utilizan los VFC del 
enfoque ráster. Finalmente, antes de las 
conclusiones, se da un ejemplo práctico de 
aplicación del enfoque ráster. 


Planteamiento matemático 


En este apartado se repasa la formulación de las 
ecuaciones de Saint-Venant bidimensionales, 
se aplica al sistema una discretización en 
volúmenes finitos en el caso general de mallas 
no estructuradas y finalmente se muestra cómo 
se simplifica la formulación cuando se adapta 
el esquema a problemas en los que es posible 
el uso de mallas estructuradas y homogéneas. 

Hay que decir que la presentación de 
este apartado se ha hecho desarrollando 
la formulación hasta el último nivel, el 
de programación, donde las expresiones 
algebraicas resultantes finales tienen gran 
sencillez. Se ha hecho de esta manera para 
mostrar la simplificación a la que se puede llegar 
en la implementación del enfoque ráster en el 
problema hidrodinámico del flujo en lámina 
libre en 2D. Se han dejado sin explicación los 
conceptos teóricos subyacentes que pueden 
encontrarse en las referencias bibliográficas 
que se dan. 


Ecuaciones de Saint-Venant 
bidimensionales 


Las ecuaciones de Saint-Venant bidimen- 
sionales, denominadas en inglés “Shallow 
Waters equations”, se pueden presentar en 
forma conservativa y notación vectorial de la 
siguiente manera: 


U+VE=H!+H 0) 


donde t es el tiempo; U, el vector de variables 
de flujo; F, el tensor de flujo; H*, la componente 
motriz del término independiente, y H”? es 
la pendiente de fricción, que responden a las 
expresiones: 


h hu hv 
y? 
U=|hu|; F= hug huv ] 
hv 2 


Q) 


donde h es el calado o profundidad del agua, u, 
la velocidad en la dirección del eje de abscisas 
y v la velocidad en la de las ordenadas; g, la 
aceleración de la gravedad; S,,, la pendiente 
del fondo en la dirección de las abscisas, Soy 
en el de las ordenadas y S,, es la pendiente de 
rozamiento en la dirección de las abscisas, y S y 
en el de las ordenadas: 


donde n es el coeficiente de rugosidad de 
Manning. La ecuación (1) representan los 
principios de conservación de la masa y de 
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conservación de la cantidad de movimiento en 
las dos direcciones del plano x/ y. 

La ecuación consta de tres términos. El 
primero representa la variación temporal local 
de las variables hidráulicas: masa y cantidad de 
movimiento; el segundo significa la variación 
espacial de los flujos de dichas cantidades, y 
el tercero (término independiente) constituye 
la ganancia O pérdida de masa y cantidad 
de movimiento por unidad de tiempo en un 
volumen diferencial que se mueve con el fluido. 
Evidentemente, la variación de masa debe ser 
nula, por no incluirse ni la infiltración ni la 
precipitación, por lo quela primera componente 
de H' es cero. La contribución exterior a la 
cantidad de movimiento, con las hipótesis 
realizadas, tiene dos razones: la variación de 
energía potencial —reflejada en la pendiente 
del fondo— y las fuerzas de fricción con el 
contorno —reflejada en la pendiente motriz—. 
Las ecuaciones de Saint Venant son un caso 
concreto de sistema de ecuaciones diferenciales 
en derivadas parciales hiperbólico, cuasi-lineal 
y con término independiente. 


Esquema numérico en volúmenes finitos 


Aplicando el Riemann solver de Roe (Roe, 1981) 
para cada lado de la celda de la figura 1, en 
un esquema en volúmenes finitos, siguiendo 
los procedimientos mostrados en Toro (2009), 
y tratando la pendiente del fondo del término 
independiente según Vázquez-Cendón (1999), 
se Obtiene el esquema numérico de primer 
orden utilizado en este trabajo. Todo este 
proceso deductivo completo puede encontrarse 
en detalle en Bladé y Gómez (2006), donde se 
llega al siguiente esquema explícito, presentado 
sintéticamente a modo de estructura recursiva: 


yr! U” At a F SS Ú 

Sd TA (U,)+F(U;))n; - Vez 
ij=l k=1 

A, +at(H) (4) 


donde: 


1 y n representan el índice denotativo 
de celda, y el de instante de tiempo 
de integración; j representa el índice 
denotativo de celda contigua a la i. 

1¡ representa el lado compartido por las dos 
celdas i y ¡ de longitud Al,. 

At es la longitud del paso de tiempo de 
integración. 

S, es la superficie en planta de la celda. 
n= (n, ny es el vector ortonormal 
apuntando hacia afuera de la celda, siendo 
n, y 1, sus componentes en las direcciones 
de los ejes coordenados. 


e,: 


1 0 1 
e =ju+cn, | ;e,=|-cn, | ;e,=|u-=cn, (5) 
v+cn), cn), 9-cn), 
Y: 
Y= ayp,+(1-signo (,))B, (6) 


Figura 1. Discretización en volúmenes finitos de un 
dominio bidimensional. En este caso se representa 
una discretización con una malla compuesta de 
cuadriláteros irregulares. Aquí, ¡ denota el elemento 
donde se aplicará la discretización y la estructura 
recursiva, S, es la superficie que encierra el volumen 
y n, es el vector ortonormal a lado correspondiente. 


Tecnología y 
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con: 


013=37 +2 -[8(1m0)n, + A(ho)n, -(un, + on,) An] 


A) = z[(a(10) - vAh)n, -(4(1u) - uAh)n, | 


P, =-0.5cAz 
B2=0 
Bi =-0.50A2), (7) 


donde A es un operador numérico que, aplicado 
sobre una variable f, vale Af= f; - f. 


M3 = 48, +0N, +C AN si y] > Ez 


y 
r => . 
dz = UN, +0N, y Ez Si |ky|<Ez 
con: €, = máx [0,2 MD; -»y] (8) 


En este esquema, las expresiones (5) a (8) 
son evaluadas sobre los lados compartidos 
mediante las siguientes interpolaciones entre 
las variables hidráulicas de celdas contiguas: 


UC; Fu; 
Y —— 
C¡+C; 


DC +DCj 
o. E; c,=Jgh,; c;=,[gh; (9) 


c¡+C; 


Adaptación a malla de volúmenes finitos 
cuadrangulares 


En este apartado se describen los cambios 
necesarios para la adaptación de la formulación 
propuesta de (4) a (9) para una malla 
estructurada y homogénea de VFC al enfoque 
ráster. 

En una malla de VFC todas sus celdas son 
cuadradas y homogéneas, y están dispuestas 
por filas y columnas, de manera que cualquier 
celda del dominio de integración tiene un 
elemento superior denotado por el subíndice 


“N”, uno inferior denotado por el subíndice 
“S”, uno a la derecha denotado por el subíndice 
“E” y uno a la izquierda denotado por el 
subíndice “0”. Los lados del volumen finito *I” 
adyacentes a cada uno de los volúmenes que 
lo rodean son subíndices en minúsculas por 
“np”, *e”, “5” y “0”, respectivamente (figura 2). 
La susodicha adaptación consiste en tomar los 
siguientes valores: 


N,=4 
0 1 
l,=Ax jm, =| [;l,=Ay;n.=| |; 
1 0 
0 -1 
l, =Ax;n, = ¡l, =Ay;n, = (10) 
- 0 


y aplicarlos a las expresiones (4) a (9). 

En la sección “Apéndice de expresiones 
algebraicas” del final de este trabajo, se dan 
todas las expresiones simplificadas resultantes 
de substituir (10) en (5) a (9). Todas ellas son de 
tipo algebraico y tienen una marcada sencillez. 
En su conjunto, sólo dependen de las cotas 
de los volúmenes finitos contiguos y de sus 
variables hidráulicas. 


=>. 


Ax 


Figura 2. Discretización en volúmenes finitos de un 
dominio bidimensional de VEC usada en el enfoque ráster. 
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Finalmente, resultará muy útil definir 
el vector de flujo mediante la siguiente 
transformación del sumatorio sobre los cuatro 
lados de (4): 


dh y 


4 k=3 
an Sej0Jorto)oy Hrs], 
k=1 


hxyuy hgUs 


hy uz, +0.52/%, - huz +0.5gh% Ax 


hyUyOy hsusos 
hgUp hi, Uy 
+ h¿UEDE - h¿ Uy Vo Ay 
hoz +0.5ghg ¿07 +0.5gh5 
1 0 1 
Yin Us +FY2n| Cn [+ Y3n Y, 
D, HC, 0 UE, 
1 0 1 
+Yis Us ul Yas Cs + Yas Ys Ax 
DC, 0 D, +C, 
1 0 1 


1 0 1 
(11) 


+1 |Mo-Co|+Y, 10 |+y uy, +Cy ||Ay 
0 0 0 


D 0 Do 


Teniendo en cuenta el vector de flujo (11) 
y los demás vectores originales, el esquema 
general de actualización temporal (4), útil para 
cualquier tipo de malla, queda adaptado para 
su utilización en el contexto del enfoque ráster 
mediante la siguiente estructura recursiva: 


h n+1 h n dh n 0 

hu| [mul LE lam at gns; | (12) 
28, 

hv ; hv ] dhv a ghSs, h 


Condiciones de contorno 


En la gran mayoría de casos, la cartografía 
utilizada para los estudios de inundabilidad y 
de clasificación en función del riesgo potencial 
de balsas presenta grandes zonas donde 
no se dispone de información altimétrica 
susceptible de ser utilizada para la creación de 
los correspondientes MDE. Los VFC o celdas 
afectados por esta falta de información se 
denominan de tipo NODATA. La razón de ello 
es el ahorro que se tiene cuando no se comparan 
zonas donde se supone, en principio, que no 
debe circular agua. 

El tratamiento que se hace sobre los 
NODATA y los contornos delimitados por las 
edificaciones es el del establecimiento de una 
condición de contorno tipo pared material, por 
donde el flujo no puede pasar. Desde el punto 
de vista numérico, cuando a un VFC se le va 
a aplicar (12) para actualizar en el tiempo las 
variables y cuando alguno de los volúmenes, 
N, E, S 0 O es de tipo NODATA, se mantiene 
el mismo esquema, pero asignando a los 
NODATA los siguientes valores: 


hi =hj 5 =-4];05 =-0) (13) 


j 
donde ahora ¡ es el subíndice de VFC tipo 
NODATA, 1, el subíndice del VFC adyacente 
al ¡ con el que comparte el lado; h, el calado; 
u, la velocidad en el eje de las abscisas, y v es 
la velocidad en el eje de las ordenadas en el 
instante n. 

La gracia de (13) está en el hecho de que se 
anulan los promedios (9) de las velocidades en 
el lado compartido, siendo nulo el flujo a través 
del mismo, pero manteniendo inalterado el 
esquema (11)-(12). 

La presencia de volúmenes NODATA adya- 
centes a un VFC da lugar a 15 combinaciones 
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posibles, que generan los correspondientes 
tipos de celda. En la figura 3 se muestra un 
esquema de todas las posibles combinaciones 
entre volúmenes normales de tipo 1 y los 
NODATA. El tipo de combinación pasa a ser 
una nueva propiedad de todas las celdas —que 
describe la estructura topológica de malla— 
y debe subministrarse con el mismo formato 
ráster de entrada de datos de las demás 
propiedades hidráulicas. 

La existencia de VFC de tipo NODATA 
da la posibilidad de utilizar esta figura para 
simular el comportamiento del flujo al entorno 
de algunas estructuras, como los pilares de los 
puentes. Ello puede llevarse a cabo siempre 
y cuando sea posible componer la sección 
del pilar mediante cuadrados NODATA de la 
precisión del MDE utilizado (longitud del lado 
de la celda). 

El tratamiento que se hace sobre los 
contornos rectilíneos de los márgenes de la 
serie cartográfica (digamos la primera y última 
filas del ráster, y primera y última columnas) 
que constituye el modelo, es tratado como 


condición de contorno de tipo salida libre de 
agua: 


h; a hy E u; = pide v; = v; (14) 


donde ahora el subíndice ¡ representa cualquier 
VFC perteneciente al contorno de la serie 
cartográfica e 1 representa el subíndice de 
cualquier VFC adyacente al contorno con el 
que comparte un lado. 


Ejemplo ilustrativo 


Como ejemplo ilustrativo se propone el 
estudio de un tramo del río Tajo (España) 
de unos 20 km, cuyo modelo digital de 
elevaciones (MDE) tiene el sistema de 
referencia horizontal UTM del uso 30-norte 
del Instituto Geográfico Nacional de España. 
En este sistema de referencia, el ámbito viene 
dado por el recuadro de X mínima: 389 522 
m; de X máxima: 397 802 m; de Y mínima: 
4 417 623 m, y de Y máxima: 4 429 603 m. Por 
lo tanto, el ancho de recuadro es de 8.280 km 


q 


Figura 3. Tipos de volúmenes con condiciones de contorno tipo Pared material en algunos de los volúmenes adyacentes. En gris 


se representa un VEC de tipo NODATA. Los 15 tipos representan todas las combinaciones posibles. La topología de la malla 


queda definida por un solo número, que pasa a ser una propiedad más de la celda. 
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y el de alto 11.980 km. Si la dimensión de la 
celda/VFC es de 5 x 5 m, entonces el número 
de columnas del ráster será de 1 656 y el de 
filas de 2 396, dando lugar a un número total 
de celdas o VFC de 3 967 776, de los cuales el 
número total de celda/VFC con información 
de cota es 717 796 y sin dicha información 
(celda/VFC de tipo NODATA) es 3 249 980. 

En el enfoque ráster, la presencia de celdas 
tipo NODATA tiene una gran influencia sobre 
el tiempo total de computación. Para evaluar 
la importancia de dicha influencia, se han 
implementado dos MDE. En el primer MDE 
—denominado Modelo con NODATA—, el 
número de VFC susceptible de ser “mojado” 
por el agua es de 717 796, el 18.1% del total 
de celdas del MDE (figura 4). En el segundo 
MDE —denominado Modelo sin NODATA—, 
dicho número pasa a ser de 3 967 776, el 100% 
de las celdas, porque se ha substituido el valor 
NODATA de la celda por una cota arbitraria 
suficientemente alta de 600 m como para no 
dejarse “mojar” por el agua (figura 4). De 
esta forma, los resultados hidráulicos son 
equivalentes. 

El tiempo que tarda una computadora en 
resolver un problema se denomina tiempo 
de computación. En el enfoque ráster, este 
tiempo se reparte entre el tiempo de escaneo 
(se entiende por escaneo el paso por todos los 
VFC en cada instante de tiempo, para ver si 
existe información de cota, es decir, se trata del 
periodo de tiempo utilizado para ver si cada 
VFC es de tipo NODATA o no; hay que notar 
que cuantificar este tiempo resulta complicado 
y que en el enfoque vectorial, este tiempo no 
existe, dado que todos los elementos, si existen, 
tienen información geográfica) y el tiempo de 
cálculo (tiempo utilizado para la actualización 
de las variables hidráulicas, según la estructura 
recursiva (12)). 

El número de celdas tipo NODATA presentes 
en un MDE tiene suma importancia, porque 
el tiempo de escaneo se alarga de manera 
enorme (más que de forma proporcional) y, en 
consecuencia, el tiempo de computación. Dado 
que el flujo de agua nunca sale del cauce con 


información de cota, el tiempo de cálculo en 
ambos modelos es exactamente el mismo. 

Por otro lado, si se usa alguna estrategia de 
paralelización en el código de programación, 
el número de celdas tipo NODATA presentes 
en el MDE no tiene por qué afectar de forma 
tan fuerte al tiempo de computación, porque 
el tiempo de avance depende del VFC, que 
requiere mayores cálculos matemáticos; es 
decir, aquel con agua que tiene los cuatro 
volúmenes del contorno lleno de agua. 

El ensayo consistió en calcular el tiempo de 
computación necesitado para la computación 
de 8 h de simulación, con la entrada de un 
caudal constante de 1 700 m*/s desde la 
sección de aguas arriba del río. El coeficiente de 
rugosidad de Manning fue establecido en 0.035 
para todos los VEC. 

Las condiciones iniciales de la sección de 
aguas arriba del tramo estudiado han sido 


Cota (m) 
MM 122172 
EA] 442.809 
EU] 453.789 
EY 464.770 
EA] 475.750 
A] 486.387 
E] 497.368 
MM 508.348 
E 519.328 


ES] NODATA 


Figura 4. Modelo digital de elevación de 20 km del río Tajo, 
España, que ha servido de base para el modelo hidráulico. 
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impuestas en régimen rápido, para poder 
asegurar la entrada de caudal propuesto. En 
concreto, la velocidad en el eje de las ordenadas 
se corresponde al doble de la celeridad en cada 
punto (figura 5). 

El resultado final del modelo hidráulico que 
se obtiene al cabo de 8 h se muestra en la figura 
6, donde puede verse que el agua no ha llegado 
hasta el final del tramo de estudio. 

Aunque los resultados del modelo hidráu- 
lico, como los mostrados en la figura 6, son 
muy interesantes, el objetivo principal de este 
trabajo es el de presentar el enfoque ráster 
y también demostrar su mayor eficiencia de 
cálculo, comparado con el enfoque vectorial. 
Por lo tanto, la comparativa en cuanto a 
tiempos de computación es el resultado más 
importante. 

La comparativa sobre tiempo de compu- 
tación entre el enfoque vectorial y el enfoque 
ráster se ha llevado a cabo resolviendo el 
mismo problema con tres programas de cálculo 
diferentes, codificados en dos lenguajes de 
programación, en dos sistemas operativos 
y en dos diferentes ordenadores. Las tres 
combinaciones de hardware/software fueron las 
siguientes: 


Figura 6. Sobre el MDE de fondo, estado de la lámina en el 


instante 28 800 s (8 h) de iniciada la simulación. 
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Figura 5. Condiciones iniciales y de contorno constantes en la sección de aguas arriba del tramo de río analizado. 
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1. 


La combinación 1 consistió en la 
computación, con el programa IBER, 
que dispone de un motor de cálculo 
escrito en código Fortran y que trabaja 
sin  paralelización. Dicho programa 
resuelve el problema del flujo en lámina 
libre en 2D, según el enfoque vectorial 
y por lo tanto está adaptado a cualquier 
tipo de malla, bien sea estructurada o no 
estructurada. Para el caso de estudio en 
cuestión, se ha implementado una malla 
estructurada, en que cada celda del MDE 
se corresponde con el consiguiente VEC. 
Así, el modelo implementado tiene tantos 
VFC como celdas tipo con información de 
cota. El caso es que el programa desconoce 
que la malla es estructurada, de manera 
que debe consultar la base de datos de la 
topología para todos los VEC a cada paso 
de tiempo de integración. Como ya se ha 
dicho anteriormente, esta estrategia hace 
lento el proceso de cálculo. El ensayo se 
ha hecho en un ordenador de 64 bits de 
última generación, que dispone de un 
CPU modelo Intel core (TM) i5 CPU 750 
2.8 GHz y de una memoria RAM de 12 GB. 
La combinación 2 consistió en la utilización 
de otro programa también escrito en código 
Fortran y sin paralelización. Este programa 
resuelve el problema del flujo en lámina 
libre en 2D, según el enfoque ráster, y 
por lo tanto, sólo resuelve los problemas 
con malla estructurada y homogénea, 
siguiendo la filosofía Do one thing and 
do it well, de acuerdo con la información 
original del MDE sin interpolación. Esta 
combinación ha sido implementada en 
el mismo ordenador en que lo ha sido la 
combinación 1. 


3. La combinación 3 consistió en el uso de 
un tercer programa escrito en código 
C-CUDA —extensión del lenguaje C para 
programación, en paralelo de la GPU de las 
placas gráficas del tipo NVIDIA— (Sanders 
y Kandrot, 2010). La  paralelización 
consistió solamente en lanzar un thread 
de actualización temporal de las variables 
hidrodinámicas para cada VFC. En este caso, 
el tiempo de escaneo resulta irrelevante, 
porque todos los thread lanzados empiezan 
por preguntar si el VFC es de tipo 
NODATA, con lo que el porcentaje de VFC 
tipo NODATA presente en un problema no 
afecta de manera proporcional al tiempo 
de escaneo. El ensayo se ha hecho en un 
ordenador diferente, también de 64 bits 
de última generación, que dispone de una 
CPU modelo Intel (R) Core (IM) 17 CPU 
950 O 3.07 GHz y de una memoria RAM 
de 11.8 GB. La placa gráfica utilizada fue la 
Geforce GTX 580 de NVIDIA. 


Las tres combinaciones de hardware / software 
han dado los resultados que se presentan en el 
cuadro 1 en cuanto a tiempo de resolución. 

A la vista de los resultados mostrados en el 
cuadro 1, caben destacar los siguientes puntos: 


e La combinación 1, trabajando con el 
modelo sin NODATA, no dispone de su 
correspondiente resultado. Ello se debe a 
que la experiencia indica que el enfoque 
vectorial con más de tres millones de 
volúmenes finitos con un programa sin 
paralelización resulta inviable desde el 
punto de vista operativo, por consumir 
demasiado tiempo de computación y, por 
consiguiente, no se ha realizado este ensayo 
numérico. 


Cuadro 1. Tiempos de computación usados en los cinco ensayos numéricos realizados. 


Modelo con NODATA (minutos) Modelo sin NODATA (minutos) 
Combinación 1 1035 - 
Combinación 2 106 402 
Combinación 3 85 158 
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Observando los resultados de las tres 
combinaciones cuando se resuelve el 
modelo con NODATA, puede verse la 
mejora sustancial que supone el enfoque 
ráster frente al enfoque vectorial, viéndose 
reducido el tiempo de computación en 9.76 
veces con la combinación 2 y en 12.18 veces 
con la combinación 3. 

Del análisis de los mismos resultados 
analizados en el punto anterior, sorprende 
la poca diferencia de tiempos existente entre 
las combinaciones 2 y 3. Ello significa que 
el tiempo de CPU requerido para el escaneo 
y Cálculo lineal de los 700 000 VFC es del 
mismo orden de magnitud que el tiempo 
de GPU requerido por el cálculo lineal del 
thread más lento (VFC con agua y con los 
cuatro del contorno también “mojados”) de 
todos los calculados en paralelo. 

La paralelización del problema planteado 
como enfoque ráster es el que presenta 
una mayor rapidez de cálculo de todas 
las posibles combinaciones y además es la 
estrategia que mejor soporta la reducción 
de VFC de tipo NODATA en el MDE, como 
puede verse al comparar los valores de la 
última fila del cuadro. 


Conclusiones 


Se enumeran las conclusiones del trabajo: 


Se estableció una clasificación de los 
problemas de flujo en lámina libre en dos 
dimensiones de acuerdo con el formato 
de la información disponible y requerida 
de los sistemas de información geográfica 
(SIG): enfoque vectorial y enfoque ráster. 

Se presentó el sistema de ecuaciones de las 
aguas poco profundas de Saint-Venant en 
dos dimensiones y se describió el método 
de resolución utilizado, que se basa en el 
esquema explícito de primer orden en 
volúmenes finitos de Godunov con el 
Riemann solver de Roe y con el tratamiento 
del término independiente propuesto 
por Vázquez-Cendón. Dado que dicho 


esquema sirve tanto para problemas con 
malla estructurada o no, el problema puede 
considerarse dentro del enfoque vectorial. 
Se establece una adaptación del método 
propuesto en el punto anterior a pro- 
blemas que se describen con una malla 
estructurada y homogénea de cuadriláteros 
que hace equivaler —directamente y sin 
interpolaciones— cada celda del MDE 
con el correspondiente volumen finito 
cuadrangular (VFC). Las dimensiones 
del VEC serán las de la celda de base del 
MDE. Por lo tanto, el algoritmo puede 
considerarse dentro del enfoque ráster. 

En el enfoque ráster, el hecho de que a cada 
instante de tiempo de integración no sea 
necesaria la consulta de la base de datos de 
la topología de la malla por ser conocida de 
antemano —precisamente porque se trata 
de una malla estructurada y homogénea—, 
hace que el proceso numérico mejore de 
forma enorme su eficiencia de cálculo 
tanto por requerir menor tiempo de 
computación como menor memoria RAM 
de almacenamiento de la geometría. 

La sencillez de planteamiento del 
enfoque ráster se mantiene también en 
el establecimiento de las condiciones 
de contorno. Se proponen dos tipos 
condiciones de contorno por defecto: 
pared material sin flujo para los VFC 
tipo NODATA —celdas sin información 
de cota— y salida libre de agua —para 
el contorno de la serie cartográfica del 
MDE—. 

El enfoque ráster reduce en gran medida 
el tiempo de computación respecto del 
enfoque vectorial. En el caso ilustrativo 
en que no se utiliza la estrategia de 
paralelización, se necesitó 9.76 veces 
menos de tiempo y cuando se utilizó 
paralelización, 12.18. 

La mejor estrategia de implementación 
de un programa de resolución del flujo 
en lámina libre en 2D en cuanto al tiempo 
de computación es el enfoque ráster con 
paralelización, porque, por un lado, puede 
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escribirse el código con cierta facilidad 
y sencillez, y por el otro, porque es la 
estrategia que mejor soporta el incremento 
de las dimensiones del problema. 


Apéndice de expresiones algebraicas 


En este apéndice se muestra la particularización 
de las expresiones (5) a (9) —del denominado 
enfoque vectorial— al caso de implementación 
del enfoque ráster, objeto del presente trabajo. 
Por particularización se entiende la substitu- 
ción y simplificación de dichas expresiones, 
tomando N, = 4 los siguientes valores: 


l, = Ax L, = Ay 
0 il 

n= n= 
1 0 
a AR l, =Ay 
0 -1 

n, = n, = 
-1 0 


Las expresiones algebraicas que se pre- 
sentan aquí están ordenadas para facilitar su 
implementación en un programa de natu- 
raleza lineal, es decir un thread. En caso de 
paralelización del problema, donde se envía 
cada VFC como thread diferente, las expresiones 
se calculan de esta forma secuencial. El orden 
de las expresiones se ha hecho mediante una 
distribución en matriz de 2 x 2, una casilla para 
cada uno de los cuatro lados del VFC. 


Variables interpoladas en los lados del VFC 


Cg+C; Co+C] 
Vectores propios 
1 0 1 
€ > Y, 1€2n =| 745 |+€3n = Us 
Vi EC, 0 V,=C 
1 0 1 
€, =|U¿+C, €, 0 ¡€32 =[U, Co 
Ve Co vo 
1 0 1 
€1,=| Us  |/€,,=|€, [35 Us 
v¿=C, 0 v¿+C, 
1 0 1 
€¡, =[U,=C, |-€2,= 0 ¡€3, =|U,¿+C, 
v, 7Co Vo 
Valores propios 
Mi» = Vr + £. Mm. => U, + Co 
Man = Vi Me > U, 
Nan = VW, > C5 Mae =4,-C, 
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Factores alfa Factores épsilon 
E = máx|0,A,, -Aj,Ay-=A 
Oy > “nt de + [Iv > hyv; =V,, (My = h; )] bl El Ñ q 1.) 
Sn €2n máx [0,2,, == MN, My = Mz» ] 
1 
Ay - [yin - hu, 4, (My E h,)] Ez, = MÁX [0,2, ET Ma] 
n 
e mot o A] €, = Máx [0,A,, — Aj, Ag —»o] 
Ñ £7, = máx[0, Az, ¿Ag =2ge] 
a tha + [go ha 0, -1)] £3, = máx [0,A3, — 21, Mg — Ago] 


e 


1 
One = — [Mv hy; -v, (Mg -h,)| 
Se 87, = máx[0,A7, - Aj,Ag=2os] 


e,,= máx[0,2;,—Aj,As- . 


hh, 1 p 
0, 30 Vte - hu, —0,(hig —h)] 87, = máx[0,Ay, — Aj, As — 25] 
e 
€, = Máx]/|0,A, —A,,A7-A 
07,9 414 2 Tigwg huy (5 ho] O 
2 2, £7, = máx[0,%7,—2y,A9=20] 
1 d 
O); = [sus hu, -u, (»; - h,)] £y, = máx [0, A, -MA9= Mao | 
s 
h¿-h, 1 
054 Us - hyv¡ vs (Mg —h)] Factores 
ho =h, 1 ' 
Oj [Moto - hu, 4, (Mo = h)| [Ay | SE n= En [y¿] Si [dy,|= 61, 
- 2e, Pin = . le = . 
1 Ejy Sl [Az] < 1, Ej, Sl |. < Elo 
a», =-—|h,vo -hyv,-v (h, —h 
Ñ al pd de ol o 5) [2.7] si [A-| > € |r2.] si |.) 22 
Ed 1 Pan = io Po = TN 
Az, = E Y +3, Poo - huy -u, (Ko -)] €2n Sl | 2n| < €2n Eze Sl | el < Ez. 
[Az] sI [Az | > €3, [3] si Mao =t30 
2 P3n = . Pz. = . 
a Factores beta Ez Sl Piar < 3, Eze Sl [Az < Eze 
E 
E =-0.5 E =-0.5c, (2, - 
P Bi, CS, (2n 21) Br Co (ze 21) p,| si Da, Sta, y, si rol = 810 
E 0 0 Vis - 
E Bar Ba. Ñ €1s si Pr, | <Exs á E10 si AA 
E =0.5c, (z,,- =0.5c, (2, 
: Ban CS, (Zn 21) Bz. Co (Ze 2,) le jet ies. [ao] Si 7, [207, 
a Pos =3 . Pao = : 
Z By, =-0.50, (2; -21) Br, = -0.5c, (Zo -2,) €25 SÍ [A),|<8z 29 Sl | “Ez 
o 
3 Ba, =0 Bao =0 [Az] si [Az¿|>83, [30] si [30] >E3, 
Z Pas = 7 30 =7 . 
Ba, =0.5c, (25 -2;) Baz, =0.5c, (Zo 21) Ez, Si |rz,|<E%, £30 Si [Az,]<E3 
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Factores gamma 


Yin = Ou Pr +(1-signo (A ))Br, 
Yan = Pan +(1- signo (A, ))Ba, 
Yan = Uan Pan +(1-signo (Ay, ))Bs, 
Yo = 0% ¿P1, +(1-signo(A.))B. 
Yo = %eP20 +(1-signo(2,.))B2. 

) 


Y3e = Az¿Pze + (1 = signo (A, )Bs. 


Ys = AP is + (1 - signo (1.))Br, 
Yos = UP; +(1 - signo (A) 
) 


Y3s = Ags¿P3s + (1 al signo (As, 


Y30 F AzoP3o + (1 = signo (Az, JB, 
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Abstract 


SOLER-GUITART, J., BLADÉ, E., BOFILL-ABELLÓ, J. +? GAMAZO, P. “Raster approach” 
to the hydrodynamic problem of 2D free surface flow. Water Technology and Sciences (in 
Spanish). Vol. IV, No. 4, September-October, 2013, pp. 77-92. 


Geographical information systems traditionally provide information using two formats—aster 
and vector. Each one of these formats is handled according to their independent algorithms. 
Consequently, any type of hydraulic study performed in a territory contains information that 
can be classified by these two formats. LIDAR is the most precise technology used to generate 
Digital Elevation Models (DEM) for large areas, with a spatial resolution of 1 to 2 m and a 
minimum precision in height. It gives the cartographic information its raster characteristic, 
with the cell as its basic unit of information. Consequently, this kind of information is very 
suitable to the generation of structured quadrilateral meshes to numerically model free surface 
flow in two dimensions. The main objective of this work was to present an adaptation of 
the classic explicit first order finite volumes scheme to a finite square volumes (FSV) mesh 
which, directly and without interpolations, makes each cell in the DEM (or raster) equal to 
the corresponding FSV. Consequently, the dimensions of the ESV are those of the base DEM 
cell. This is presented herein as a raster approach to the hydrodynamic problem of free surface 
flow in two dimensions. With the raster approach, the database for the mesh topology does not 
need to be consulted for each instant of time and, therefore, the efficiency of the numeric process 
is greatly improved. Consequently, larger-sized problems can be addressed. An additional 
objective was to compare both approaches by solving an illustrative example. 


Keywords: two-dimensional modelling, shallow water, parallelization, finite volumes. 
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